Probing a liquid to glass transition in equilibrium 
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We use computer simulations to investigate the static properties of a simple glass-forming fluid in 
which the positions of a finite fraction of the particles has been frozen in. By probing the equilibrium 
distribution of the overlap between independent configurations of the liquid, we find strong evidence 
that this random pinning induces a glass transition. At low temperatures, our numerical findings are 
consistent with the existence of a random first order phase transition rounded by finite size effects. 



Experiments allow to measure the increase of the vis- 
cosity of supercooled liquids approaching the glass transi- 
tion over more than 15 orders of magnitude. Despite this 
large range, there is at present no consensus whether this 
phenomenon is controlled by an underlying phase tran- 
sition, or whether relaxation times progressively increase 
down to zero temperature [1]. There is not even agree- 
ment on the microscopic mechanisms at work, despite 
the fact that glassy dynamics is commonly observed in 
a large variety of materials, from simple liquids, network 
forming liquids, to soft and biological matter [2j. Because 
supercooled liquids inevitably fall out of equilibrium at 
the experimental glass temperature T g , amorphous solids 
form without encountering any singularity. 

The existence of a thermodynamic transition at a fi- 
nite temperature T K < T g is highly debated \2\. Such 
a singularity was discussed long ago by Kauzmann [3] 
who pointed out that experimental data for the config- 
urational entropy extrapolate to zero at Tk , the "Kauz- 
mann temperature" . Physically, this would imply that 
for T > Tk there exist exponentially many (in the 
number of particles) equilibrium states (neglecting vi- 
brations), whereas below Tk this number becomes sub- 
exponential and the system forms an "ideal" glass. This 
view is supported by the analysis of the temperature de- 
pendence of the relaxation dynamics, which seems to be 
consistent with the existence of a temperature T > 
at which the relaxation time diverges [J]. This makes 
sense if thermodynamic and dynamical singularities oc- 
cur at the same temperature, To sa Tk, but so far To 
and Tk have been determined only by empirical fitting 
formula and uncontrolled extrapolations of the behavior 
of macroscopic observables. These procedures evidently 
leave much room for debates on the existence of such a 
transition [51 [3]. 

The existence of a true phase transition is appealing 
from a theoretical perspective, since a number of analyt- 
ical results obtained within various approximate schemes 
or idealized limits suggest the existence of a "random 
first order transition" (RFOT) occurring in glass-forming 
materials [3 [5]. This transition results from the tem- 
perature evolution of a complex free energy landscape 
characterized by a large number of metastable states. It 
provides an exact realization of Kauzmann's entropy cri- 
sis. However, the mean-field character of the RFOT, the 
technical difficulties posed by finite-dimensional fluctu- 
ations 9], and the indirect connections between theory 



and experiments explain why RFOT is only one among 
several theoretical scenarios for the glass transition [2J. 

In the following we demonstrate that it is indeed pos- 
sible to perform a detailed equilibrium study of a liquid- 
glass transition whose nature is equivalent to the one 
possibly occurring at Tk in bulk glass-formers [THl E] • 
This allows us to probe numerically, for the first time, 
the nature of the phase transition and the microscopic 
properties of the glass phase at thermal equilibrium. In 
particular, we show that glass formation can be located 
without using uncontrolled extrapolations, and that equi- 
librium glasses can be produced. 

We consider a 50:50 binary mixture of harmonic 
spheres of diameter ratio 1.4 at constant density 

p = 0.675, which we study using molecular dynamics (see 
SM for methodological details). This quasi-hard sphere 
system has an onset temperature around T on ps 10, and a 
mode-coupling temperature T mct ~ 5.2 [T3]. All numbers 
are expressed in appropriate reduced units (see SM). To 
efficiently sample the equilibrium thermodynamic prop- 
erties of the system, we use replica exchange molecular 
dynamics [14] . We carefully ensure thermalization and 
equilibrium sampling by requiring that all particles move 
several times across the entire simulation box for each 
state point and that all replicas properly explore the con- 
figuration space. 

Our central idea is to induce a glass transition by in- 
creasing the strength of a random pinning field in a ther- 
malized dense liquid [Till [UJ. We first equilibrate the 
liquid at a given temperature T, before freezing perma- 
nently the position of a certain set of particles. Among 
various other pinning geometries [TTJ Q21 [TB] , here we ran- 
domly select a finite concentration of particles, c. Work- 
ing at constant total number density p, we then study the 
equilibrium properties of the model in the (c, T) phase di- 
agram. For each state point, we average the results over 
independent realizations of the random pinning field |15j . 
Intuitively, pinned particles constrain the available phase 
space of the remaining fluid particles, thus impacting 
their static and dynamic properties. Indeed, we find that 
the dynamics slows down dramatically with c at constant 
T, see Fig. SMI, in agreement with earlier studies jTT|[T7T- 

EH]. 

Freezing particles at the positions they occupy at equi- 
librium is a key ingredient to this study, since this does 
not perturb the interaction part of the Hamiltonian, but 
only constrains the available phase space. Indeed, any 
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FIG. 1: (a) and (b): Large spheres represent pinned particles 
(rescaled in size by a factor 0.5), small dots are the superpo- 
sition of the position of fluid particles obtained from a large 
number of independent equilibrium configurations at T = 4.8 
for c = 0.0625 (a) and c = 0.1875 (b). The total number 
of particles is N = 64. Panel (c): Average overlap (q) as a 
function of concentration c of pinned particles for N = 64 and 
different temperatures. 



ensemble averaged correlation function measured in the 
pinned liquid takes the same value as in the bulk sys- 
tem at c = P]3 HO] . Theoretical studies have recently 
revealed two additional consequences [TUl [21] . First, in- 
creasing c affects the configurational space in the same 
way as decreasing temperature does for the bulk liq- 
uid [TO]- In particular, the RFOT occurring at for 
c = becomes a transition line Tk(c) > with equiv- 
alent properties. Second, the pinning field only induces 
a glass transition if temperature is low enough so that 
metastable states and glassy dynamics are well-developed 
at c = 0. The glass transition line Tk(c) ends at a sec- 
ond order critical point [10], and no transition exists for 
temperatures above this point. Note that other types of 
pinning fields 22-24J produce qualitatively different re- 
sults, because they perturb the Hamiltonian and possibly 
affect the nature of the glass transition |21) . 

In Fig. [TJ we provide a qualitative illustration of our 
strategy. We show how an equilibrium glass state is 
obtained at constant temperature, T — 4.8, if the con- 
centration c of pinned particles is increased. In these 
snapshots the large spheres are the pinned particles in 



a typical realization of the disorder whereas the small 
black dots stem from the superposition of a large num- 
ber of independent, equilibrium configurations visited by 
the fluid particles. For c — 0, the dots look like mist ho- 
mogeneously filling the simulation box. For c = 0.0625, 
panel (a), the fluid particles are not much constrained 
and still have access to a large number of distinct config- 
urations. Thus, the dots form fuzzy clouds. An increase 
to c = 0.1875, panel (b), condenses the dots into well- 
defined patches, which represent the highly constrained 
positions occupied by the free particles. This qualita- 
tive observation illustrates that increasing c leads to a 
collective localization of the fluid particles. We empha- 
size that during the simulations, all fluid particles diffuse 
and explore the entire simulation box. However, particles 
move in such a way that collective density fluctuations 
are frozen. Thus, the system is in a glass state character- 
ized by a frozen amorphous density profile [2] , but which 
can nevertheless be studied in equilibrium conditions be- 
cause our replica exchange simulation algorithm permits 
single particle diffusion. 

We now study in more detail how the system makes 
the transition between the fluid and glass states shown 
m Fig.fl] The images suggest that the number of avail- 
able states, and thus the configurational entropy, consid- 
erably decreases with c. A quantitative determination of 
the configurational entropy is however difficult and has 
several shortcomings [2] • Therefore we use a microscopic 
order parameter to characterize the transition from fuzzy 
clouds to small patches seen in Fig. [T] An appropri- 
ate quantity is the overlap q a p measuring the degree of 
similarity between two arbitrary configurations a and ft, 
which has been used in spin glass models displaying an 
RFOT [7 . In practice, we discretize space into small cu- 
bic boxes of linear size 0.55, and define = 1 if box i 

in configuration a is occupied by a particle, and rif 1 ' = 
if not. Then, 

, N b 
k i—l 

where the sum runs over the Nb boxes which do not con- 
tain pinned particles. By definition q aa = 1, while a 
small overlap is obtained for independent configurations 
(t7 ran d ~ 0.11 on average for an ensemble of independent 
configurations). In the remainder of the paper, we char- 
acterize the transition through a detailed analysis of the 
statistical properties of the overlap for a broad range of 
control parameters. 

In Fig. Ilk, we show the c-dependence of the average 
overlap, (q) = {q a p), where the brackets stand for ther- 
mal and disorder averages. Above the onset temperature, 
(g) increases gradually with c. For T < 8.0, the growth 
remains modest at low c, but this initial regime is fol- 
lowed by a rapid increase in the range (q) s» 0.25 — 0.4. 
Finally, if T is decreased even further this rapid growth 
occurs for lower values of c and becomes sharper. At 
T — 4.8 it is sufficient to reach c « 0.11 to abruptly lo- 
calize the fluid particles, (q) > 0.5, while such a large 
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f} P{q)qdq at high T shown in Fig. lip. A qualita- 
• different behavior is observed at lowT, see Fig. [2b. 





FIG. 2: Probability distribution of the overlap, Eq. B, for 
different values of the concentration of pinned particles at 
high and low T (panels (a) and (b) respectively). Note the 
presence of a double peak structure at low temperatures. 



overlap requires three times as many pinned particles, 
c « 0.35, near the onset temperature. 

These observations show that at low temperatures 
the glass state is reached at a sharply defined c-value, 
whereas pinning acts quite smoothly at high T. The In- 
dependence of (q) is consistent with the emergence of 
a discontinuous jump at low temperatures, T < 8.0, but 
the relatively small system size shown in Fig.[lp, N = 64, 
will obviously smear out such a discontinuity. To over- 
come this difficulty we have determined the fluctuations 
of the overlap and in Fig. [2] we show the overlap distri- 
bution function, 



P(.q) = (8(q-q af t)). 



(2) 



At high temperatures, T = 13 in Fig. [2po, P(q) evolves 
smoothly from distributions peaked at small q for small 
c to distributions peaked at large q values at large c. 
This is consistent with P(q) becoming a delta-function 
for N —¥ oo for all c, and with the smooth increase of 



(?) 

tively i 

A narrow peak is still present for small and large c, but 
P(q) is bimodal for intermediate c. The presence of two 
peaks implies that it is equally probable that two inde- 
pendent thcrmalizcd configurations are either very sim- 
ilar or very different. The former situation is favored 
at larger c because too few distinct configurations exist, 
while the latter holds at small c when the pinning field 
is not strong enough to prevent the system from explor- 
ing a large configuration space. Bimodal distributions of 
the order parameter can be interpreted as the phase co- 
existence between the low-g liquid and the high-g glass 
phases, and the two-peak structure of P(q) is suggestive 
of a first-order transition for the order parameter (g), 
rounded by finite size effects. 

We emphasize that the above conclusions do not re- 
sult from any kind of extrapolation, but follow directly 
from equilibrium measurements of a microscopic order 
parameter on both sides of the transition. To estimate 
more precisely the location of this putative transition, 
we have measured various moments of the distribution 
P(q). From the second moment (q 2 ) = J Q P(q)q 2 dq, we 
define the static susceptibility, 



X (C,T)=N(1-C) [(q 2 )-(q)' 



(3) 



The curves presented in Fig. [3^ show that at high tem- 
peratures x nas a mild c-dependence, whereas it devel- 
ops a well-defined peak for T < 7.0 whose height and 
location respectively increases and shifts to smaller c if 
T decreases. The observed growth of the peak x(c, T) at 
low T is a direct evidence of increasing static correlations 
in the bulk system [TTJ [55] . The location of this peak al- 
lows to estimate the value of the critical concentration at 
the given temperature. As demonstrated in Fig. |3]d, an 
even more accurate location is obtained by considering 
the skewness j(c, T) of the distribution P(q), 

((g-(g)) 3 ) 
{(q-(q)?f 2 ' 



(4) 



since 7 crosses zero when the distribution is symmetric, 
i.e. precisely at coexistence. From the location of the 
peak in x(c, T) and the zero-crossing point of 7(0, T), we 
can locate the fluid-glass coexistence line, as shown in 
Fig.[3p. These two quantitative estimates become reliable 
if T < 8.0, i.e. when bimodal distributions are present, 
and we see that they agree with each other within er- 
ror bars (estimated using the jack-knife method). In the 
context of Fig. [3p, the existence of a Kauzmann temper- 
ature at c = would rely on extrapolating our finite c 
transition points (determined without extrapolation) to 
the limit c — > 0. 

Altogether, our results show that a glass transition is 
induced by random pinning, with nontrivial thermody- 
namic signatures occurring at a sharp value of the ran- 
dom pinning field that are consistent with the existence of 
an equilibrium random first-order transition [10J. How- 
ever, since our results have been obtained for a rather 
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FIG. 3: (a) The static susceptibility x, Eq. pi, as a function of 
concentration c for different temperatures, (b) c-dependence 
of the skewness 7 for different T. (c) The equilibrium phase 
diagram determined from the position of the peak in the sus- 
ceptibility x an d from the root of the skewness 7 showing 
the location of the fluid and glass phases. Also included are 
curves of constant relaxation time as determined from the 
self-intermediate scattering function. 



FIG. 4: System size dependence of the overlap distribution 
function for T — 5.5 and three values of c below (green), 
near (red), and above (blue) the transition for N = 64 (lines) 
and -/V = 128 (lines with symbols). The peak susceptibility 
increases from x « 1-8 for N = 64 to x ~ 2.7 for N = 128 
(data not shown). 



small system size, TV = 64, they do not establish its ex- 
istence in the thermodynamic limit. Unfortunately, the 
parallel tempering algorithm we use becomes less effi- 
cient when N is increased. While we cannot access very 
low temperatures for larger N, we can, however, investi- 
gate the iV-dependence in a limited range. Our results 
for larger systems confirm that no sharp transition exists 
for T > 8.0, while we obtained clear signs of enhanced 
bimodality and static susceptibility when T < 6.0, as 
shown in Fig. [4] for T — 5.5. The small dip in P(q) for 
N = 64 and intermediate c values becomes more pro- 
nounced for N = 128, while the peaks at low and large 
overlaps become sharper. Such system-size dependence 
is again typical of first-order transitions. A more exten- 
sive finite-size scaling analysis would be useful, but it is 
at present beyond the computational means. 

In contrast to conventional studies of the glass transi- 
tion, our approach allows us to study both liquid and 
glass phases at thermal equilibrium, it does not re- 
quire uncontrolled extrapolations, and we characterize 
the glass formation by means of a microscopic order pa- 
rameter. Our work thus demonstrates the feasibility of 
systematic equilibrium studies of the nature of the liquid- 
to-glass transition. Our approach is very general and 
can easily be applied to any type of glass-former. Also, 
by using optical tweezers it will be possible to use this 
method to study colloidal glasses and thus to produce ex- 
perimentally what is often considered as "impossible-to- 
reach ideal glass states" [26] . Therefore our work opens 
the door for a new generation of direct, systematic stud- 
ies of the nature of the glass state and for the production 
of novel amorphous materials. 
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Supplementary Material 



I. MODEL AND DETAILS OF THE 
SIMULATIONS 

The model we have been using is a 50:50 binary mix- 
ture of N harmonic spheres pQ. All particles have the 
same mass m and particles i and j separated by a dis- 
tance Tij have an interaction given by 



V(r ij ) = -(l-(r ij /a ij ) 2 ) 



(5) 



Here an = 1.0, a i2 = 1.2, and a 2 2 = 1-4. We will use an 
as the unit of distance, \J man/e as the unit of time and 
10~ 4 e as unit of energy, setting the Boltzmann constant 
ks = 1-0. The density we have been using is 0.6749715. 

Before we discuss the details of the procedure how we 
have pinned the particles, we recall that pinning parti- 
cles in a liquid configuration that has been equilibrated 
at a temperature T pin does not perturb the equilibrium 
properties of the remaining liquid at the temperature 
Tiiq = Tpi n if one averages over a sufficiently large num- 



ber of realisations of pinned configurations [2r|4] . There- 
fore it is not necessary to re-equilibrate the pinned sys- 
tem after the pinning. If instead T p ; n is different from 
the simulation temperature Ti; q of the liquid the resulting 
Kauzmann line Tr-(c, Tn q ) will depend on T p ; n as well [5]. 

To prepare the system at different concentrations c of 
pinned particles we have proceeded as follows. First we 
have used a standard molecular dynamics simulations to 
equilibrate a liquid (in the bulk) at the target tempera- 
ture Tii q . For this we have integrated the equations of mo- 
tion using the velocity form of the Verlet algorithm with 
a time step of 0.02. Let us call one of the so obtained bulk 
configurations Cbuik- Subsequently we have to determine 
which ones of the particles have to be pinned. Although 
chosing these pinned particles randomly from the Cbuik 
configuration is certainly a possibility, such a choice can 
lead to systems in which many pinned particles are clus- 
tered or that have relatively large regions in which there 
are no pinned particles. These fluctuations make that 
the transition at the Kauzmann point becomes smeared 
out at intermediate temperatures and hence makes it dif- 
ficult to locate this transition accurately. Therefore we 
have chosen a scheme that makes that the typical dis- 
tance between pinned particles, and therefore the size of 
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the remaining liquid regions, is relatively uniform. In 
practice we have generated an equilibrated liquid config- 
uration of c • N particles, at density p — 0.672, using 
a molecular dynamics simulation with periodic bound- 
ary conditions. Let us call a typical configuration of this 
simulation C pin . Since this is a liquid like configuration, 
all particles have a similar number of nearest neighbors 
and also a well defined nearest neighbor distance from 
each other. Subsequently we have scaled this configura- 
tion such that the box size becomes the same as the one 
of configuration Cbuik- Then we have searched for each 
particles of the C P i n configuration the closest particle in 
the Cbuik configuration. These closest particles have then 
been pinned permanently. In short this procedure allows 
to select from the Cbuik configuration c • N particles that 
can be pinned permanently and that have a relatively 
unifom distance from each other. 

Typically we have averaged our results over around 
5 independent realisations of the pinning field if c was 
small and over up to 20 realisations when c was large, i.e. 
around the Kauzmann point. The total amount of CPU 
time used for obtaining the presented results is about 
900,000 hours. 

II. RELAXATION TIMES 

In order to determine the relaxation time of the pinned 
system as a function of temperature and concentration c 
we have used standard molecular dynamics simulations 
with a system of N = 1000 particles (step size 0.05). 
From these runs we have determined the self intermediate 
scattering function F s (q,t), where q is the wave- vector. 
Using for q the value of 5.5, which corresponds to the po- 
sition of the maximum in the static structure factor, we 
have defined the relaxation time r via F s (q, r) = e _1 . In 
Fig. [5] we show the temperature dependence of t(T, c) in 
an Arrhenius plot. From this graph one recognizes that 
for c = 0, i.e. the bulk system, the data curves slightly 
upwards, i.e. the system is a somewhat fragile glass- 
former, as discussed in Ref . [1] . An increase of c has two 
effects: 1) The dynamics is slowing down very strongly 
(by several orders of magnitude at intermediate and low 
T) and this slowing down becomes more pronounced if 
T is decreased. 2) At the highest c considered the re- 
laxation dynamics becomes Arrhenius-like, i.e. the data 
falls on a straight line. This observation is in qualitative 
agreement with the theoretical prediction by RFOT [TU] 
which shows that this aspect of the theory seems to be 
compatible with the simulation data. (We mention that 
a similar decrease of fragility with increasing c has been 
observed in Ref. [7].) 

Since Fig. [5] demonstrates that the relaxation time is 
a strong function of c, one can expect that the relevant 
temperatures of the glass-former, such as the onset tem- 
perature To, the mode-coupling temperature T c , and the 
Kauzmann temperature Tk will also depend on c and 
will increase if c is increased. As a consequence it should 
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FIG. 5: Arrhenius plot of the relaxation time for different 
value of the concentration c of pinned particles. 



become possible to cross the line Tk{c) by doing equilib- 
rium simulations at low T by increasing c [8] and in 
the main text we show that this crossing, and hence the 
investigation of the occurring transition, can indeed be 
done in equilibrium. 

Since the standard molecular dynamics slows down 
quickly upon approaching the Tr-(c) line, we have used 
the parallel tempering algorithm [91 [TU] using 24 repli- 
cas. The smallest difference in temperature between two 
neighboring replicas was 0.3, which has allowed for a 
good overlap of the two potential energy distributions. 
Attempts to switch two neighboring replicas have been 
made every 50,000 time steps. We have checked that each 
parallel tempering run has indeed reached equilibrium by 
following any given replica and making sure that all tem- 
peratures have been sampled sufficiently. A typical path 
of a replica in temperature space is shown in Fig. [6] for 
our lowest studied temperature. 

The parallel tempering algorithm does indeed allow 
to equilibrate the system down to T = 4.8 even if the 
concentration c of pinned particles is large. This can be 
recognized from Fig. [7] where we show the mean squared 
displacement (MSD) of the particles (distinguishing small 
and large particles) as a function of time. To calculate the 
quantity we have followed a given replica in temperature 
space and considered only time intervals at which this 
replica was at the target temperature Tn q . The figure 
shows that at sufficiently long times the MSD becomes 
very large, thus indicating that the particles do indeed 
move through the box also at high values of c, even if 
their relative arrangment does not change much, i.e. one 
is in the glass state. Or stated otherwise: Each particle 
explores all the small patches shown in Fig. la of the 
main text. 
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time step 



FIG. 6: Time dependence of a typical trajectory of a replica in 
temperature space. The temperatures are labeled 0-23, with 
being the target temperature Tn q . The example shown is 
for Tii q = 4.8 and the concentration of pinned particles is 
c = 0.14. 
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FIG. 7: Mean squared displacement of the particles as a func- 
tion of time using the parallel tempering dynamics. The circles 
and crosses correspond the the large and small particles, re- 
spectively. T = 4.8 and c = 0.14. 
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